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Algorithmic aspects of multicanonical simulations 
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Monte Carlo (MC) simulations of many systems, in particular those with conflicting constraints, can be con- 
siderably speeded up by using multicanonical or related methods. Some of these approaches sample with a-priori 
unknown weight factors. After introducing the concept, I shall focus on two aspects: (i) Opinions about the 
optimal choice of weight factors, (ii) Methods to get weight factor estimates, with emphasize on a multicanonical 
recursion. 



1. INTRODUCTION 

One of the questions which ought to be ad- 
dressed before performing large scale computer 
simulations is "What are suitable weight factors 
for the problem at hand?" . It has been expert 
wisdom for quite a while 0, and became widely 
recognized in recent years [2-8], that MC simu- 
lations with a-priori unknown weight factors are 
feasible and deserve to be considered. 

The design of suitable weights requires to un- 
derstand or guess relevant physical features of the 
system(s) under (numerical) investigation. This 
is already exemplified by the work of Metropolis 
et al.[^. To sample the important configurations 
of a canonical ensemble at temperature T = (3^^ 
efficiently, the Boltzmann weight is recommended 



(1) 



where Ea is the energy of configuration a. 
Namely, it is straightforward to show that this 
weight minimizes the sampling error for the 
canonical (normalized) energy density P{E): 



(AP(_E))^ = minimum. 



(2) 



However, even when studying canonical statisti- 
cal mechanics the Boltzmann weight is not a sat- 
isfactory recipe for all cases. It can be false that 
the weight factor (1) samples the important con- 
figurations. A well-studied case is the sampling of 
configurations with interfaces Quite gen- 

erally, instead of (2) the sampling error for the 



actually calculated observable needs to be opti- 
mized. 

A second reason, called "dynamic" in re- 
lies on the fact that the error estimate (2) holds 
for independent configurations. To sample with a 
weight like (1), one has to rely on a Markov chain. 
This introduces an autocorrelation time which de- 
pends on the weights chosen. In this sense the ar- 
gument leading from (1) to (2) is "static" Q. It 
can happen that statically most efficient weights 
are rendered inefficient by the dynamics, whereas 
statically worse weights can yield a favorable dy- 
namics. For instance, in systems with conflict- 
ing constraints encouraging experience has been 
made with weights and expanded ensembles 
Ik -13 1 which allow to "refresh" the system in the 
disordered phase. It seems worthwhile to study 
this concept also for less difficult situations, like 
second order phase transitions. 

In the next section I consider the question of 
optimal weights. The third (and major) section 
focuses on how to get working estimates of the 
weights in the first place. Brief conclusions are 
given in section 4. 



2. WHICH WEIGHTS? 

First, I consider the static aspects. For the rel- 
ative sampling error of the spectral density n(E) 
we get 



'y^^{An{E) / n{E))'^ = minimum 



(3) 
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by using 

w{a) = l/n{Ea) . (4) 

The result is an uniform energy probability den- 
sity. For a large number of applications this has 
proven to be a robust choice, for instance when 
dealing with first order phase transitions Pp^. 
Instead of the energy, weights in other thermody- 
namic variables have also been used successfully, 
for example w{a) — f{ma) 0], where rua is the 
magnetization of configuration a, or wla) = f{ba) 
Q , where b the bound variable for cluster updat- 
ing. 

In some applications one is tempted to intro- 
duce weight factors in two variables, like w{a) = 
f{Ea,ma), which would allow to re-construct 
canonical expectation values over both, a macro- 
scopic temperature and magnetic field range. 
Only for exceptional cases, when one is interested 
in very small volumes, this seems to be practi- 
cal. In most cases the interest lies in a finite size 
scaling (FSS) analysis of the infinite volume limit 
V oo. CPU time slows down worse than V"^ for 
the discussed methods. One cannot afford to re- 
place V hyV"^ , i.e. a performance worse than V^. 
In addition, a RAM increase with V'^ is unpleas- 
ant too. The art seems to be to find one combi- 
nation of variables, such that suitable weighting 
with it leads to significant improvements of the 
numerical performance. 

With respect to groundstate searches it is 
claimed in ref.||] that the weight 

w{a) = l/ £ n{E') (5) 

gives the best worst case performance in terms 
of ergodicity and pertinence. The authors define 
pertinence by the static argument that VNs in- 
dependent samples with the weight (5) are better 
than Na independent rival samples. Ergodicity 
relates to the dynamics of the Markov chain. The 
latter is system dependent and there is presently 
little understanding of the processes driving its 
slowing down. Hence, it seems to me that their 
statement about ergodicity is unwarranted. 

The idea of expanded ensembles ||l^,|ll|, par- 
ticularly attractive in their version of parallel 



tempering ||lj,|l^, tries to overcome dynamical 
slowing down by introducing moves into addi- 
tional dimensions. Whereas new weights may al- 
ter or by-pass barriers, expanded (canonical) en- 
sembles can only by-pass them, because Boltz- 
mann weights are used for the normal updates. 
A, yet untested, possibility would be to apply the 
parallel tempering idea to unconventional ensem- 
bles. This may allow to combine the advantages 
of both approaches. 

3. RECURSIVE WEIGHT ESTIMATES 

Weight factors like (4) or (5) are a-priori un- 
known. This introduces some additional compli- 
cations, unknown to canonical simulations. Be- 
fore the simulation can be performed, a working 
estimate of the weight factor has to be known, i. e. 
an approximation for the practical purpose of per- 
forming the MC simulation. Deviations by factors 
around ten may well be acceptable, because the 
accuracy of the simulation is hardly affected by 
this. The errors of the simulation remains in any 
case statistical, because the actually used weight 
factors are exactly known. On the other hand, 
deviations from the desired weights by factors of, 
say, order 10^*^ (well the order of improvements) 
are certainly no longer acceptable. Such large dis- 
crepancies would entirely obstruct the basic ideas. 

In some cases, for instance first order phase 
transitions ||l6|, FSS allows reasonable work- 
ing estimates on the basis of already simulated 
smaller systems. Also constraint simulations have 
successfully been used |]l7| , |l8| . For other situa- 
tions, like spin glasses y|, this does not work. In 
the remainder of this section I will first explain 
the difficulties of the most straightforward recur- 
sion. Subsequently, I give an efficient solution, 
which essentially is a simplified and improved pre- 
sentation of an result of rcf. |l9) . 

Desired is a working estimate of the weight fac- 
tor (4). Others can be treated similarly. It is 
advisable to start any recursion in the disordered 
phase of a system, where the system moves most 
freely under MC updates. Therefore, let us as- 
sume for the zeroth simulation. 

w'^(a) = 1 for all a. (6) 
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The most obvious recursion goes now as fol- 
lows: Simulation n, {n — 0, 1, 2, ...) is carried out 
with the estimate w"{a) and yields the histogram 
H'^{E). Estimate n + 1 for the weight factors is 
then given by 



w"(a) 



(7) 



This recursion has a number of difficulties 

(i) What to do with histogram entries H^\E) — 

or small? 

(ii) Each H"{E) calculation starts with zero 
statistics. Assume, someone has given us 
the exact result and we use it for w^{a): 
The next estimate w^{a) will be worse. A 
noisy left-over of w'^{a). 

(iii) Initial weights w{a) — 1 correspond to tem- 
perature infinity and are bad in the limit 
Ea —>■ En-iin, approached towards zero tem- 
perature. 



To avoid (i), one may replace 
H{E) H{E) = max [/iq, H{E)] 



(8) 



where /iq is a number < 1. The recursion derived 
in the following overcomes (ii). 

of 



Let us first discuss the relationship |20 



ic 



the weight factors with microcanonical tempera- 
ture (3{E), because it turn out to be advantageous 
to derive the final recursion in this quantity. It 
holds 

«;"+i(a) = e-^(^») = g-/3(E„)iJ„+o(Eo) (9) 

where S{E) is the microcanonical entropy and, 
by definition. 



p{E) 



dS{E) 
dE 



This determines the fugacity function a{E) up to 
an (irrelevant) additive constant. Consider, for 
example, the case of a discrete minimal energy e. 
We may choose 



and the identity S{E) = I3{E) E - a{E) imphes 
S[E) - S{E - e) = P{E)E - l3{E - e){E - e) - 
a{E) +a{E-e). Inserting e (3{E - e) = SiE) - 
S{E~e) yields 

a{E ^ e) ^ a{E) + [(3{E - e) - (3{E)] E (11) 

and a{E) is fixed by defining Q!(i?inax) = 0. In 
summary, once P{E) is given, a{E) follows for 
free. The starting condition ^ becomes (other 
j3^{E) choices are of course possible) 



j3''{E)=cP{E)=Q. 



(12) 



With H defined by (g) we translate now eqn.(7) 
into an equation for (3{E). Subscripts o are used 
to indicate that those quantities are not yet the 
final estimators from the n*'' simulation. Let 



w"{E) 
H^{E) 



where the (otherwise irrelevant) constant c is in- 
troduced to ensure that Sq^^{E) can be an esti- 



mator of the microcanonical entropy. It follows 
Si;+\E) = -\nc + S'^iE) +\nH"-lE). Inserting 
this relation into {nw gives 



P^+\E)^P-{E) + 

[\nH''{E + e) - \nH"{E)]/e 

As estimator of the variance follows 



(13) 



+ 



m{E + €) ' w{E) ' 

where c' is an unknown constant and we have 
re-introduced the original histograms to empha- 
size (and use) that the variance is infinite when 
there is zero statistics. The statistical weight for 
Pq'^^{E) is inversely proportional to its variance. 
Choosing a convenient constant, we get 

H'^iE + e) H"'{E) 



9o(.E) = 



W^iE- 



W{E) 



(14) 



Note that gl}{E) = for H''{E + e) = oi 
H'''^{E) — 0. The n*'* simulation was carried out 
using P^{E). It is now straightforward to com- 
bine f3o~^^{E) and P'^iE) according to their re- 
spective statistical weights into the desired esti- 
mator: 



(3{E) = [S{E + e) - S{E)] /e 



(10) 



p-+\E)^g-{E)p-{E)+g^iE)/3^+'iE), (15) 
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where the normaUzed weights 
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9oiE) 



and g'\E) = 1 - gl\E) 



are determined by the recursion 

g^+\E)^g-{E)+g-^{E), g\E) = 0. (16) 



We can ehminate (3q^^{E) from eqn.(|l5|) by in- 
serting its definition ( p^ ) and get 
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(17) 



Notice that the sole purpose of using H, instead of 
H has become to avoid ln(0). The weight g^iE) 
of such contributions is zero. The major advan- 
tage of (p7|), compared with (|^), is that it keeps 
already assembled statistics. This solves (ii) and 
eases (iii). Frequent iterations are allowed, imply- 
ing increased stability and decreased CPU time 
consumption. In addition one may implement a 
suitable 13{E) guess for E fi'min- 

Finally, eqn.(^^ can be converted into a re- 
cursion for ratios of weight factor neighbors. We 
define 



i?"(£;) 



w"(£;) 



w"(£; + e) 
and get the recursion 



(18) 



R^+'{E) = R'\E) 



H''{E + e) 
m{E) 



ao(E) 



(19) 



4. CONCLUSIONS 



Sampling of broad energy and other unconven- 
tional distributions is technically feasible. Stati- 
cal slowing down can be overcome, as is well es- 
tablished for first order phase transitions. Dy- 
namical slowing down is more difficult, but fur- 
ther algorithmic improvements (find the relevant 
combination of variables!) appear to be possible. 
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